#################################################################################
###Do File: Climate Adaptation and Conflict Mitigation:The Case of South Sudan###
#################################################################################
library(doBy)
library(countrycode)
library(foreign)
library(dplyr)
library(lfe)
library(sandwich)
library(DataCombine)
library(lmtest)
library(stargazer)
library(interplot)
library(plm)
#Set to relevant working directory
setwd("~/OneDrive - Indiana University/FromGoogle/SouthSudanProject/Code_12_20_24/")
#Read in the dataset
env.dat <- read.csv("ssdat92024.csv")
summary(env.dat)


############################
############################
###Table 4.2 (OLS Models)###
############################
############################

#####################
###Baseline Models###
#####################
##Civil War
lm.1.cw <- lm(acled_all_cw~lagclim_env_fs+lagconf_pb+laglogIDPs+
                laglogged.pop+laglogp_avg+lagt_c_avg+
                lagacled_all_cw+t+factor(month),
              data=env.dat)
summary(lm.1.cw) 
lm.1.cw.r <- coeftest(lm.1.cw, vcov. = vcovCL(lm.1.cw, cluster = env.dat$gid, type = "HC0"))
lm.1.cw.r

##Social Conflict
lm.1.sc <- lm(acled_all_sc~lagclim_env_fs+lagconf_pb+laglogIDPs+
                laglogged.pop+laglogp_avg+lagt_c_avg+
                lagacled_all_sc+t+factor(month),
              data=env.dat)
summary(lm.1.sc)
lm.1.sc.r <- coeftest(lm.1.sc, vcov. = vcovCL(lm.1.sc, cluster = env.dat$gid, type = "HC0"))
lm.1.sc.r



#################################
###General Preparedness Models###
#################################
##Civil War
lm.2.cw <- lm(acled_all_cw~lagclim_env_fs.gp+lagconf_pb+laglogIDPs+
                laglogged.pop+laglogp_avg+lagt_c_avg+
                lagacled_all_cw+t+factor(month),
              data=env.dat)
summary(lm.2.cw) 
lm.2.cw.r <- coeftest(lm.2.cw, vcov. = vcovCL(lm.2.cw, cluster = env.dat$gid, type = "HC0"))
lm.2.cw.r

##Social Conflict
lm.2.sc <- lm(acled_all_sc~lagclim_env_fs.gp+lagconf_pb+laglogIDPs+
                laglogged.pop+laglogp_avg+lagt_c_avg+
                lagacled_all_sc+t+factor(month),
              data=env.dat)
summary(lm.2.sc)
lm.2.sc.r <- coeftest(lm.2.sc, vcov. = vcovCL(lm.2.sc, cluster = env.dat$gid, type = "HC0"))
lm.2.sc.r

#Results hold for when we use the NA dataset (env.dat) POTENTIALLY REPORT


###############################
###Community Building Models###
###############################
##Civil War
lm.3.cw <- lm(acled_all_cw~lagclim_env_fs.cb+lagconf_pb+laglogIDPs+
                laglogged.pop+laglogp_avg+lagt_c_avg+
                lagacled_all_cw+t+factor(month),
              data=env.dat)
summary(lm.3.cw) 
lm.3.cw.r <- coeftest(lm.3.cw, vcov. = vcovCL(lm.3.cw, cluster = env.dat$gid, type = "HC0"))
lm.3.cw.r

##Social Conflict
lm.3.sc <- lm(acled_all_sc~lagclim_env_fs.cb+lagconf_pb+laglogIDPs+
                laglogged.pop+laglogp_avg+lagt_c_avg+
                lagacled_all_sc+t+factor(month),
              data=env.dat)
summary(lm.3.sc)
lm.3.sc.r <- coeftest(lm.3.sc, vcov. = vcovCL(lm.3.sc, cluster = env.dat$gid, type = "HC0"))
lm.3.sc.r


##Export to Latex
#RSE
stargazer(lm.1.cw.r,lm.2.cw.r,lm.3.cw.r,lm.1.sc.r,lm.2.sc.r,lm.3.sc.r)
#Regular
stargazer(lm.1.cw,lm.2.cw,lm.3.cw,lm.1.sc,lm.2.sc,lm.3.sc)

############################
############################
###Table 4.3 (GMM Models)###
############################
############################

#Create binary month indicators
env.dat$month1 <- ifelse(env.dat$month==1,1,0)
env.dat$month2 <- ifelse(env.dat$month==2,1,0)
env.dat$month3 <- ifelse(env.dat$month==3,1,0)
env.dat$month4 <- ifelse(env.dat$month==4,1,0)
env.dat$month5 <- ifelse(env.dat$month==5,1,0)
env.dat$month6 <- ifelse(env.dat$month==6,1,0)
env.dat$month7 <- ifelse(env.dat$month==7,1,0)
env.dat$month8 <- ifelse(env.dat$month==8,1,0)
env.dat$month9 <- ifelse(env.dat$month==9,1,0)
env.dat$month10 <- ifelse(env.dat$month==10,1,0)
env.dat$month11 <- ifelse(env.dat$month==11,1,0)
env.dat$month12 <- ifelse(env.dat$month==12,1,0)

#Convert data to pgmm formal
m.pgmm <- pdata.frame(env.dat, index=c("gid", "t"))

#####################
###Baseline Models###
#####################
##Civil war
gmm.1.cw <- pgmm(acled_all_cw~lagclim_env_fs+lagconf_pb+laglogIDPs+
                   laglogged.pop+laglogp_avg+lagt_c_avg+lagacled_all_cw+month2+                     
                   month3+month4+month5+month6+month7+month8+month9+month10+month11+month12| 
                   lag(acled_all_cw, 2:7),
                 data = m.pgmm, effect = "individual", transformation= "ld", model = "twostep")
summary(gmm.1.cw)

##Social conflict
gmm.1.sc <- pgmm(acled_all_sc~lagclim_env_fs+lagconf_pb+laglogIDPs+
                   laglogged.pop+laglogp_avg+lagt_c_avg+lagacled_all_sc+month2+                     
                   month3+month4+month5+month6+month7+month8+month9+month10+month11+month12| 
                   lag(acled_all_sc, 2:7),
                 data = m.pgmm, effect = "individual", transformation= "ld", model = "twostep")
summary(gmm.1.sc)

#################################
###General Preparedness Models###
#################################
##Civil war
gmm.2.cw <- pgmm(acled_all_cw~lagclim_env_fs.gp+lagconf_pb+laglogIDPs+
                   laglogged.pop+laglogp_avg+lagt_c_avg+lagacled_all_cw+month2+                     
                   month3+month4+month5+month6+month7+month8+month9+month10+month11+month12| 
                   lag(acled_all_cw, 2:7),
                 data = m.pgmm, effect = "individual", transformation= "ld", model = "twostep")
summary(gmm.2.cw)

##Social conflict
gmm.2.sc <- pgmm(acled_all_sc~lagclim_env_fs.gp+lagconf_pb+laglogIDPs+
                   laglogged.pop+laglogp_avg+lagt_c_avg+lagacled_all_sc+month2+                     
                   month3+month4+month5+month6+month7+month8+month9+month10+month11+month12| 
                   lag(acled_all_sc, 2:7),
                 data = m.pgmm, effect = "individual", transformation= "ld", model = "twostep")
summary(gmm.2.sc)

###############################
###Community Building Models###
###############################
##Civil war
gmm.3.cw <- pgmm(acled_all_cw~lagclim_env_fs.cb+lagconf_pb+laglogIDPs+
                   laglogged.pop+laglogp_avg+lagt_c_avg+lagacled_all_cw+month2+                     
                   month3+month4+month5+month6+month7+month8+month9+month10+month11+month12| 
                   lag(acled_all_cw, 2:7),
                 data = m.pgmm, effect = "individual", transformation= "ld", model = "twostep")
summary(gmm.3.cw)

##Social conflict
gmm.3.sc <- pgmm(acled_all_sc~lagclim_env_fs.cb+lagconf_pb+laglogIDPs+
                   laglogged.pop+laglogp_avg+lagt_c_avg+lagacled_all_sc+month2+                     
                   month3+month4+month5+month6+month7+month8+month9+month10+month11+month12| 
                   lag(acled_all_sc, 2:7),
                 data = m.pgmm, effect = "individual", transformation= "ld", model = "twostep")
summary(gmm.3.sc)

##Export to Latex
stargazer(gmm.1.cw,gmm.2.cw,gmm.3.cw,gmm.1.sc,gmm.2.sc,gmm.3.sc)


####################################
####################################
###Table 4.4 (Interaction Models)###
####################################
####################################

#####################
###Baseline Models###
#####################
##Civil War
lm.1.cw.d <- felm(acled_all_cw~lagclim_env_fs*disaster+lagconf_pb+laglogIDPs+
                    laglogged.pop+laglogp_avg+lagt_c_avg+
                    lagacled_all_cw+t|month|0|gid,
                  data=env.dat)
summary(lm.1.cw.d) 
##Social Conflict
lm.1.sc.d <- felm(acled_all_sc~lagclim_env_fs*disaster+lagconf_pb+laglogIDPs+
                    laglogged.pop+laglogp_avg+lagt_c_avg+
                    lagacled_all_sc+t|month|0|gid,
                  data=env.dat)
summary(lm.1.sc.d)

#################################
###General Preparedness Models###
#################################
##Civil War
lm.2.cw.d <- felm(acled_all_cw~lagclim_env_fs.gp*disaster+lagconf_pb+laglogIDPs+
                    laglogged.pop+laglogp_avg+lagt_c_avg+
                    lagacled_all_cw+t|month|0|gid,
                  data=env.dat)
summary(lm.2.cw.d) 
##Social Conflict
lm.2.sc.d <- felm(acled_all_sc~lagclim_env_fs.gp*disaster+lagconf_pb+laglogIDPs+
                    laglogged.pop+laglogp_avg+lagt_c_avg+
                    lagacled_all_sc+t|month|0|gid,
                  data=env.dat)
summary(lm.2.sc.d)

###############################
###Community Building Models###
###############################
##Civil War
lm.3.cw.d <- felm(acled_all_cw~lagclim_env_fs.cb*disaster+lagconf_pb+laglogIDPs+
                    laglogged.pop+laglogp_avg+lagt_c_avg+
                    lagacled_all_cw+t|month|0|gid,
                  data=env.dat)
summary(lm.3.cw.d) 

##Social Conflict
lm.3.sc.d <- felm(acled_all_sc~lagclim_env_fs.cb*disaster+lagconf_pb+laglogIDPs+
                    laglogged.pop+laglogp_avg+lagt_c_avg+
                    lagacled_all_sc+t|month|0|gid,
                  data=env.dat)
summary(lm.3.sc.d)

##Export to Latex
#RSE
stargazer(lm.1.cw.d,lm.2.cw.d,lm.3.cw.d,lm.1.sc.d,lm.2.sc.d,lm.3.sc.d)


####################################
####################################
###Figure 4.3 (Interaction Plots)###
####################################
####################################

#############################
###Create FELM ME Function###
#############################
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.05, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### defining function to get marginal effects for specific levels of the treatment variable
  getmfx_high_low <- function(betas, data, treatment, moderator, treatment_val) {
    betas[treatment] * treatment_val + betas[paste0(treatment, ":", moderator)] * data[, moderator] * treatment_val
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting  marginal effects for high and low cases of treatment variable
  data[, "marginal_effects_treatment_low"] <- getmfx_high_low(coef(regression_model), data, treatment, moderator, quantile(data[,treatment], 0.05))
  data[, "marginal_effects_treatment_high"] <- getmfx_high_low(coef(regression_model), data, treatment, moderator, quantile(data[,treatment], 0.95))
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  # Marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_point(aes(y = (marginal_effects) )) + 
    geom_errorbar(aes(ymin = (ci_lower), ymax = (ci_upper ) ), width = 0.2) +  
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    scale_x_continuous(breaks = c(0, 1), labels = c("Regular Climate Stress", "High Climate Stress")) +  # Customize x-axis labels
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of", treatment_translation, "\non", dependent_variable_translation)) #+  # Update y-axis label
  # ggtitle("Average normalized marginal effects")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}


#########################
###Subset NA Omit Data###
#########################
env.dat.omit <- na.omit(env.dat[,c("acled_all_cw","acled_all_sc","lagclim_env_fs","lagclim_env_fs.gp","lagclim_env_fs.cb","disaster","lagconf_pb","laglogIDPs","laglogged.pop","laglogp_avg","lagt_c_avg","lagacled_all_cw","lagacled_all_sc","t","month","gid")])
#####################################
###Calculate for All CAFSI###
#####################################
###Estimate model with NA omit data
lm.1.sc.d <- felm(acled_all_sc~lagclim_env_fs*disaster+lagconf_pb+laglogIDPs+
                    laglogged.pop+laglogp_avg+lagt_c_avg+
                    lagacled_all_sc+t|month|0|gid,
                  data=env.dat.omit)

###Create marginal effects plot
jpeg("Figure43a.jpeg", width = 6, height = 6, units = 'in', res = 500)
felm_marginal_effects(regression_model = lm.1.sc.d, 
                      data = env.dat.omit, 
                      treatment = "lagclim_env_fs", moderator = "disaster", 
                      treatment_translation = "CAFSI (t-1)", 
                      moderator_translation = "Climate Stress Level", 
                      dependent_variable_translation = "Social Conflict")
dev.off()


#############################################
###Calculate for General Prepardness CAFSI###
#############################################
###Estimate model with NA omit data
lm.2.sc.d <- felm(acled_all_sc~lagclim_env_fs.gp*disaster+lagconf_pb+laglogIDPs+
                    laglogged.pop+laglogp_avg+lagt_c_avg+
                    lagacled_all_sc+t|month|0|gid,
                  data=env.dat.omit)

###Create marginal effects plot
jpeg("Figure43b.jpeg", width = 6, height = 6, units = 'in', res = 500)
felm_marginal_effects(regression_model = lm.2.sc.d, 
                      data = env.dat.omit, 
                      treatment = "lagclim_env_fs.gp", moderator = "disaster", 
                      treatment_translation = "CAFSI with GP (t-1)", 
                      moderator_translation = "Climate Stress Level", 
                      dependent_variable_translation = "Social Conflict")
dev.off()


############################################
###Calculate for Community Building CAFSI###
############################################
###Estimate model with NA omit data
lm.2.sc.d <- felm(acled_all_sc~lagclim_env_fs.cb*disaster+lagconf_pb+laglogIDPs+
                    laglogged.pop+laglogp_avg+lagt_c_avg+
                    lagacled_all_sc+t|month|0|gid,
                  data=env.dat.omit)

###Create marginal effects plot
jpeg("Figure43c.jpeg", width = 6, height = 6, units = 'in', res = 500)
felm_marginal_effects(regression_model = lm.2.sc.d, 
                      data = env.dat.omit, 
                      treatment = "lagclim_env_fs.cb", moderator = "disaster", 
                      treatment_translation = "CAFSI with CB (t-1)", 
                      moderator_translation = "Climate Stress Level", 
                      dependent_variable_translation = "Social Conflict")
dev.off()